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Предлагается итерационный метод решения дифференциальных уравнений, имеющих дивергентный вид и описывающих од- 
номерное стационарное течение с переходом через скорость звука. Метод основан на использовании априорной информации 
о монотонном возрастании числа Маха вдоль сопла. Сравнение с другими методами проводилось на точном решении, а также 
при расчете двухфазного течения. 


Введение 

В [ 1 ] был предложен и далее в [2] уточнен итера- 
ционный метод решения одномерных стационар- 
ных уравнений газовой динамики, который за счет 
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использования априорной информации о моно- 
тонном изменении числа Маха внутри рассматри- 
ваемой области по скорости сходимости на поря- 
док превосходит метод установления на основе схе- 
мы Маккормака [3] и имеет простой алгоритм. В 


Естественные науки 


этом методе в качестве переменных используется 
полная энтальпия и энтропийная функция, а нали- 
чие в уравнениях правых частей позволяет учесть 
неравновесные процессы. 

Однако, как отмечено в [4], методы, использую- 
щие уравнения, записанные в виде законов сохра- 
нения массы, импульса и энергии, имеют преиму- 
щества перед недивергентной формой записи. 
Кроме того, запись в дивергентной форме имеет 
более простой вид по сравнению с записью, ис- 
пользуемой в [1]. Поэтому в данной статье рассма- 
тривается дальнейшее развитие метода [1]. Предла- 
гаемый ниже подход, обладая простотой метода [1] 
и несколько уступая ему в скорости, использует ди- 
вергентную форму записи уравнений. 


Математическая постановка задачи 


Рассмотрим течение в сопле с переходом через 
скорость звука. Ставится задача нахождения параме- 
тров течения р, II, Р, удовлетворяющих одномерным 
стационарным уравнениям для идеального совер- 
шенного газа при наличии неравновесных процес- 
сов. Уравнения записаны в дивергентном виде [5]: 


ААріІ _ д с іА(рѴ 2 + Р) 
Ах ’ Ах 


А'Р + АС „ 


сІАрІЖ _ АГ , 
Ах 


( 1 ) 


где: Н=Ру/(р{у—\))+1Р/7 - полная энтальпия; р, 
Л, Р, у — плотность, скорость, давление и показа- 
тель адиабаты газа; А - площадь поперечного сече- 
ния сопла; х - продольная координата, принадле- 
жащая рассматриваемой области [х а \х ь \, С) и С 2 - 
известные, в общем случае нелинейно зависящие 
от параметров газа правые части уравнений движе- 
ния и энергии, связанные с неравновесными про- 
цессами. Штрих обозначает производную по «л». 
Все величины безразмерны. Переход к безразмер- 
ным переменным производился путем отнесения: 
х — к половине размера минимального сечения, 
А — к площади в минимальном сечении тіп(И), II- 
к критической скорости II, , р - к критической 
плотности р., Р— к произведению р, III. Н- к ква- 
драту критической скорости К 2 . 

Полагается, что задана площадь сопла А(х); на 
входе в сопло заданы граничные условия в виде 
Н=Н(х а ), 5=5(х а ), где 8=Р/р у — энтропийная функ- 
ция; С 1 (х а )=С 2 (х а )= 0. В качестве априорной исполь- 
зуется информация о том, что существует только 
одна внутренняя точка х,, в которой число Маха 
М= 1; внутри рассматриваемой области число Маха 
монотонно возрастает от дозвуковой до сверхзву- 
ковой величины и не равно нулю: 

М(х)*0,хе[х а ;х ь \, 

М(х а )< 1, М(х ь ) >1, М(х,)= 1, х, е(х я ;х А ). (2) 

Вместо первого уравнения в ( 1 ) воспользуемся 
его интегралом АрІІ=0, где С - неизвестная кон- 
станта. Введем обозначения: 


К, = А(рІІ 2 + Р), К 2 = АрІІН, (3) 


( у *.) 2 

2Р 2 (у 2 - 1 ) 


( 4 ) 


При этом второе и третье уравнения в (1) при- 
мут вид: 


^ = А'Р + АС 1 , ^ = АС 
ах ах 


( 5 ) 


Обратный переход от К ь К 2 , С к переменным р, 
И ’ Р производится по следующим формулам: 


_ уК 1 + ]( уД, У 2 (у- 
С(у+1)-^{С(у + 1)] С(у + 1) 


р = (6) 

где знак +(-) соответствует сверхзвуковому (дозву- 
ковому) течению. 

Из условий (2) следует, что существует един- 
ственная точка х=х., в которой 

(М 2 -I ) 2 =тіп. 


Очевидно, что при этом выполняются соотно- 
шения: 

С = Цх,), х, = аг§(тіп {Цх),х е(х я ;х 4 )}), (7) 


АЬ 

Ах 


= 0 . 


( 8 ) 


Таким образом, если в рассматриваемой обла- 
сти известно распределение /?,(х) и К 2 (х), то соот- 
ношения (4), (7) и (6) позволяют определить кри- 
тическое сечение х,, расход С и величины р, II, Р. 

Поскольку распределение Р (х) и К 2 (х) находит- 
ся из решения задачи Коши для обыкновенных 
дифференциальных уравнений (5), то рассмотрим 
нахождение граничных условий для этих уравне- 
ний. Предположим, что расход С известен. Так как 
имеются Н(х а ), 5(х а ), С,(х а )= 0, С 2 (х а )= 0, то в сече- 
нии х а можно найти II, р и Р из соотношений: 

(^Г‘ 8 (х а ) + Цр- - н (х а )Ц у -'= О, 

Р = 4т’ Р = Р Г 3(*аУ ( 9 ) 


Подставляя найденные значения в (3), получим 
граничные условия К[(х а ) и К 2 (х а ) для уравнений (5). 


Численный алгоритм 

С учетом выше изложенного для решения си- 
стемы уравнений ( 1 ) предлагается следующий ал- 
горитм: 

1. Вводится расчетная сетка х, и сеточные функ- 
ции р„ Ц, Р„ К и , К 2І , /=0,1, ...,&, где к- число то- 
чек сетки. Обозначим через верхний индекс у 
номер итерации. 

2. Задается у— 1 , О и начальное приближение р„ Ц, 
Р„ удовлетворяющее (2). 
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3. С использованием (3) находится тах{Л 1; , 
і=0,...,к} и присваивается всем 

4. Переход к следующей итерации: у'=/+ 1 . По зна- 
чениям р„ Ц, Р і вычисляются правые части для 
уравнений (5). 

5. Из (9) находятся граничные условия Кі(х а ) и 
К 2 (х а ) для уравнений (5). 

6. Из задачи Коши для обыкновенных дифферен- 
циальных уравнений (5) находятся К и и К 2І ; ис- 
пользуется схема Эйлера второго порядка точ- 
ности. 

7. Для ускорения сходимости с использованием К и 
выполняется нижняя релаксация: 
Р$= К^+оА Я,— Щ"' 1 ), где ] - номер текущей 
итерации, К^ 1) - значение с предыдущей итера- 
ции, - будет использоваться на текущей ите- 
рации, К и - вычислено на текущей итерации, 
со - параметр релаксации, 0<п)<1. 

8. Из (7) определяется х. Так как при использова- 
нии сеточных функций значение х, будет при- 
надлежать дискретному множеству {х 0 ,...,х А .}, то 
более точное х. находится следующим образом. 
Пусть Т,=тт{/=1,...Д}. Параболическая интер- 
поляция функции Ь по трем точкам и (8) дает: 

х _ ( х і ~ х і + 1 )Тг-і + (х 2+1 -х ІЛ )Ц + (^. і — -у )Ц + \ 
7((Х/ — х /+1 )Т 71 +(х /+ | — х 7 (х 7 . [ — х І )^ 7+ і) 

Далее, используя х,, уточняем О: 

0 = 1 (х,-х І )(х,-х І+1 ) | 

{х І _ 1 -х І )(х І _,-х І+1 ) 

+1 (х,-х І _ 1 ){х,-х І+1 ) | ь (х -Х 7 .,)(х -X,) 

1 {х І -х ІЛ ){х І -х І+1 ) ,+1 (х 7+1 -х 7 . 1 )(^ +1 -^ )' 

9. Из уравнений (6) вычисляются р„ II., Р;, если 
х,<х,, выбирается дозвуковое решение, в про- 
тивном случае - сверхзвуковое. 

10. При необходимости следующей итерации про- 
изводится переход на пункт 4. 


Сравнение с точным решением 


Здесь и в следующем разделе рассматривалось 
течение в радиусно коническом сопле Лаваля, ко- 
торое с учетом обезразмеривания описывалось за- 
висимостью Д(х)=у 2 (х), где: 


3,125, х<х„ х, =-0,5-3, 25 8 іп(0, 785398) 


у(х) = 


2,125 + ^Дщх-ТхЦц хе(х 1 ;х 2 ], х 2 =х 1 +8Іп(0,785398) 

1.625 - х + 2х ъ , хе(л: 2 ;х 3 ], х 3 = -0, 625 8Іп(0, 785398) 

1. 625- ^0, 625 : -л: 2 , іеЦ,;і 4 ], х, =0,625 8Іп(0, 26 1799) 
1, 625 - 0, 625 со 5(0,26 1 799 ) + (.г - )1§(0,26 1 799), х>х 4 


х„=— 4, х 4 =2, у= 1,4, число точек сетки &=40, началь- 
ное поле параметров (пункт 2) находилось с ис- 
пользованием соотношения из [6]: 


шіі і(А) 


= М 


у + 1 
2 + (7-1)М 2 


Л 


(/+ 1)/(2( Г -1)) 


Для получения точного решения в качестве правых 
частей использовались соотношения из [7]: 


_Нр ( дг А' 8' V Р г $' 

1 Н 0 {Ы А (1 -у)8) 1-у’ 

_ 1Шр ( М > А' 8' Л 
2 Н 0 [и А (1 -у)8)’ 


здесь И=А а [Ь 2 {(х-х,)/(х а -х,)) 2 +1-Ь 2 ]\ ЗЧЭД^х-х^+І]; 
Н=Щ[{5,/т/АУ'\ А=А(х а ); /> 2 =0,5; (>,=0,2; 
Я 0 =(7+1)/(2(7-1)); 51,= 1/у. Критическое сечение 
задано: х,=1. 


Использование указанных правых частей по- 
зволяет найти точное решение в виде: 

( 

М\ 


7+1 


2 + (7-1)М- 


тіп(2Ѵ) 


N 


II = М 


Г 2(7 ~1)Я Ѵ/2 
2 + (7 -1)М 2 


Р = 


тіп(ТѴ) 
АІ 7 


Р=8рР 


На рис. 1 изображено положение точки х в процес- 
се решения: 1 - предлагаемым методом, при ш=0,1; 
2 - итерационным методом [1]; 3 — методом устано- 
вления с использованием явной схемы Маккормака 
[3]. Для всех методов решение сходится к точному ре- 
шению х.=1. Предлагаемый метод, немного уступая 
итерационному методу [1], на порядок превосходит 
метод установления по скорости сходимости. 



Рис. 1. Точное решение. Сходимость разных методов 


На рис. 2 показано влияние параметра со на схо- 
димость итерационного процесса в предложенном 
методе. Необходимо отметить, что при малых пра- 
вых частях уравнений (1) можно использовать боль- 
шое значение параметра со, однако при этом увели- 
чивается вероятность появления осцилляций. Для 
ускорения сходимости можно первоначальное зна- 
чение со изменять через несколько итераций. На 
рис. 2 номеру 1 соответствует ш=0,06; 2 - о>=0,08; 


3 — й)=0, 1 ; 4 — © = 



/ <10 

. В последнем слу- 

І >ю 


чае видно заметное ускорение сходимости. 

На рис. 3 приведен профиль сопла и получен- 
ное распределение числа Маха, которое для всех 
рассматриваемых методов сошлось к точному ре- 
шению, при этом х,= 1 . 
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Рис. 2. Точное решение. Влияние со на сходимость предло- 
женного метода 



Рис. 3. Профиль сопла Лаваля и распределение числа Маха 
вдоль него 

Сравнение для двухфазного течения 

С использованием упомянутых методов были 
проведены расчеты двухфазного неравновесного те- 
чения в сопле Лаваля. Начальное поле параметров 
(пункт 2) находилось аналогично предыдущему раз- 
делу. Конденсация, испарение, коагуляция и дро- 
бление частиц второй фазы не рассматривались. Ко- 
эффициенты межфазного взаимодействия и осо- 
бенности расчета параметров второй фазы приведе- 
ны соответственно в [1] и [8]. Использовались сле- 
дующие параметры: давление торможения 50- ІО 5 Па; 
температура торможения Г 0 =3000 К; динамическая 
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вязкость газа 5Т0“ 5 Па-с при Т<;. весовая доля второй 
фазы 0,4; число Прандтля 0,7; теплоемкость веще- 
ства второй фазы 1420 Дж/(кгК); молекулярная 
масса смеси 30 кг/кмоль; показатель адиабаты газа 
1,1; вторая фаза монодисперсная, состоящая из 
А1,0 3 ; диаметр частиц второй фазы ІО -5 м. 

Кроме этого в описанном выше алгоритме из- 
менены пункты 3 и 7: 

- пункт 3 - по формулам (3) вычисляются Щ ) и 

- пункт 7 - выполнялись релаксации: 

К$=КГ+(0(К ІГ КП И Д^ДГ+й>(Я 2 ,-ЯГ), а 
со изменялось следующим образом: 

/0,9; ] <5 
со = < 

[0,45; у >5 

Полученные результаты продемонстрировали 
сходимость к решению х,«0, 137, аналогичную пред- 
ставленной на рис. 1 . Процесс сходимости показан 
в таблице, где для разных номеров итерации у при- 
водятся разности х, на двух соседних итерациях. 


Таблица. Сходимость разных методов. Значения аЬз( хЛ 11 ~х> ф ) 


Номер итерации 

Предлагаемый 

метод 

Метод [1] 

Метод устано- 
вления [3] 

5 

3,2-10~ 2 

3,3-Ю" 2 

1,4-ТО" 2 

10 

1,3-10~ 5 

5,7-10" 5 

8,1 -10 -3 

15 

7,2-10^ 

8,5-10 7 

5,9-10~ 3 

20 

3,3-10~ 9 

1,1-Ю" 9 

2,9-10 3 

500 

2,6-10~ 9 

0 

5,7-10 5 

920 

1,2-10"'° 

0 

9 , 4 - 10 -* 


Видно, что если в методе установлении с ис- 
пользованием явной схемы Маккормака [3] до сов- 
падения 6 цифр после запятой требовалось около 
900 итераций, то в предлагаемом методе и в методе 
[1] было достаточно 15 итераций. 

Заключение 

Изложенный итерационный метод решения 
уравнений газовой динамики унаследовал простоту 
и скорость сходимости ранее предложенного мето- 
да, но в отличие от последнего использует уравне- 
ния в дивергентной форме. Это позволяет реко- 
мендовать приводимый метод для расчета одно- 
мерных стационарных неравновесных газодинами- 
ческих течений, имеющих единственный переход 
через скорость звука с монотонным возрастанием 
числа Маха вдоль рассматриваемой области. 

5. Годунов С. К., Забродин А.В., Иванов М.Я., Крайко А.Н., Про- 
копов Г.П. Численное решение многомерных задач газовой ди- 
намики. - М.: Наука, 1976. - 400 с. 

6. Лойцянский Л. Г. Механика жидкости и газа. - М.: Наука, 1987. 
- 840 с. 

7. Галкин В.М. Пример точного решении и тестовые расчеты для 
одномерных стационарных уравнений газовой динамики //Ма- 
тематическое моделирование. - 2005. - Т. 17. - № 1. - С. 3-9. 

8. Глазунов А.А., Рычков А.Д. Исследование двухфазных нерав- 
новесных течений в осесимметричных соплах // Известия АН 
СССР. Механика жидкости и газа. - 1977. - № 6. - С. 86-91. 


21 





